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re-heating scenarios for the Quark Gluon Plasma (QGP) within this approach. 

PACS numbers: 47.75.+f, 47.10.-g, 24.10.Nz, 25.75.-q 



I. INTRODUCTION 



The fluid dynamical description of the evolution of strongly interacting matter created in heavy-ion collisions, 
initially proposed and applied to describe p+p collision at low energies, was pioneered by Landau Ever since then, 
it has been successfully used to model different colliding heavy ions at a wide range of energies. Nowadays, one of the 
most intriguing and important experimental discoveries at the Relativistic Heavy Ion Collider (RHIC) in Brookhaven, 
US, the measurement of collective flow in non-central Au-f-Au collisions, demonstrates the predicting power of the 
fluid dynamical approach. 

Experimental evidence in single-particle transverse momentum distributions, like radial flow, and in the coefficients 
of the asymmetric azimuthal distribution around the beam axis, the directed transverse flow vi, the elliptic flow V2 
and the anti-flow W3, shows that thepredictions of perfect fluid dynamical models assuming initial conditions from 
the Color Glass Condensate (CGC) overestimate certain data Q. In particular, the elliptic flow, V2, surmised to 
be created in the early stage of the collision signaling an early onset of thermalization, could only be reproduced - 
using perfect fluid dynamical calculations with an initial condition of Glauber type - with thermalization time tq < 1 
fm/c, up to transverse momenta p± < 1.5 GeV 

Remarkably, the perfect fluid dynamical calculations with the CGC initial state using a realistic description of 
the dissipative hadronic corona, could still not reproduce the elliptic flow data Q. This instigates that additional 
dissipation must happen in the fluid dynamical stage: the matter created in high energy heavy ion collisions can not 
be completely described by perfect fluid dynamics with zero viscosities and without heat conduction. 

Triggered by these developments there is an increasing interest in relativistic dissipative fluid dynamics. There are 
several recent suggestions and modifications and renewed discussions of the old enigmas of relativistic viscous fluids 
@i S 9 ■ The most intensely investigated problems are the stability of the homogeneous equilibrium solutions and the 
causal propagation of perturbations. These issues are related to each other. 

There are several investigations connecting causality and stability in dynamical systems described by hyperbolic 
partial differential equations. According to these results the mathematical properties of symmetric hyperbolic equa- 
tions ensure that the propagation speed of perturbations is finite. Due to certain additional restrictions on the material 
properties, these characteristic speeds are less than the speed of light. For the so called divergence type theories ^ it 
is straightforward to determine the conditions of causality for the full set of nonlinear evolution equations. Moreover, 
for these theories causality implies stability. In case of Israel-Stewart fluids it has been proven that the symmetric 
hyperbolicit y o f the perturbation equations is equivalent to the conditions of linear stability of the homogeneous 
equilibrium [lol. [ll| . However, the original non-perturbed equations in these theories are not known to be symmetric 
(let alone causal) for arbitrary fluid states. 

On the other hand, stability also has certain implications on causality. If the homogeneous equilibrium is asymp- 
totically stable, then the causality region of the theory shall be reduced due to damping. The causality region can 
be restricted by physical characteristic speeds [H, In this sense relativistic parabolic theories are viable, provided 
one can prove the stability of the homogeneous equilibrium. Several authors argue that hyperbolic extensions of the 
Navier-Stokes-Fourier system would not have experimental co nseq uences; the essential part of dissipative relativis- 
tic hydrodynamics is the parabolic Navier-Stokes-Fourier core [IJ]. According to these arguments the proof of the 
stability is a most fundamental issue in dissipative theories. 

As it is well known, the so called flrst order theories are unstable . The stability of the homogeneous equilibrium 
in several recent, second order theories either was not investigated [a, 01) or the obtained stability conditions seem to 
be too restrictive Q. It is important to note, that the stability conditions in the Israel-Stewart theory are complicated, 
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and cannot be conceived in a natural way, see for example Eq. (70) in Ref. fTo*! . This circumstance is in strong contrast 
to the nonrelativistic case, where the thermodynamic (equivalently hydrodynamic) stability and the positivity of the 
viscosities and the heat conduction coefficient ensure the linear stability of the homogeneous equilibrium without any 
further elaboration. 

We have recently analyzed [l^ the physical reasons of instabilities in the first order theory of Eckart J^. Our 
conclusion was that stability of the homogeneous equilibrium can be restored independently of the chosen frame 
(Eckart or Landau-Lifshitz), by exploring the physical difference between momentum density and energy flux in 
the local rest frame. We have given a minimal, stable extension of the first order Eckart theory by correcting the 
traditional treatment of the energy density. Unlike in previous relativistic theories of dissipative fluids, where the 
local rest frame energy is considered as the inter nal energy, we have suggested to apply the absolute value of the 
local rest frame energy vector, e — ^ UaT^T^u'^ ■ We proved that in this case the positivity of the coefficients in 
the classical linear response theory (heat conduction coefficient and viscosities) and the conditions of thermodynamic 
stability are sufficient to avoid generic instabilities. In this theory no further conditions are needed, in full analogy to 
the non-relativistic Navier-Stokes- Fourier equations. 

In this paper our aim is to extend this approach to a set of hyperbolic equations. Our starting point is the 
fiducial equation of state between the entropy density and the energy density. In this relation we use Lorentz scalar 
combinations of the energy-momentum tensor at zero pressure and the flow four-velocity. The analysis will be carried 
out without specifying the flow frame. We derive the relaxation equations of heat flux and viscous pressures from the 
corresponding entropy production. Then we investigate a certain generalization of the Bjorken flow and show that 
the stationary solution implies vanishing heat flux. Finally we solve the subsequent equations and investigate the 
correspondence of expansion, cooling and re-heating in our approach. 



II. THERMODYNAMICS 



In this section we derive the equations of viscous fluid dynamics using the convention where the upper indices denote 
contravariant while the lower indices covariant four vectors. The metric tensor is given as Qij = diag(l, — 1, — 1, — 1) 
and all indexes fc, .. run over 0, 1, 2, 3. We use natural units, h — k = c ~ 1, except the final section. 

The projection of the energy-momentum tensor, = UiT^^ can be interpreted as the energy flux in the local 
comoving frame, the scalar projection, e = UiT^^Uj = E^Uj, comprises the local, relativistic energy density. We have 
studied the modified equation of state, s(e, n), with = EiE'' in Ref. [13] and found that it leads to a stabilization 
of the known generic instabilities. The change of variables of the entropy density was supported by arguments from 
modern non-equilibrium thermodynamics based on the Liu procedure applied to first order weakly nonlocal state 
spaces ■ The original Eckart theory considers both time-space non diagonal components of the energy- momentum 
tensor, the energy flux and the momentum density, as dissipative contributions. We have argued that using e as 
internal energy density, results in a distinction in their physical role, and restores the stability of the homogeneous 
solutions. 

Now we extend this approach by utilizing other Lorentz scalars in combinations transverse to . With the energy 
like Lorentz scalars, e^, EiE^ and TijT^^ , we construct a fiducial energy density expression containing all the dissipative 
and heat conducting terms and derive relaxation equations for them. The energy-momentum tensor can be generally 
decomposed [Ifi] into reversible and irreversible parts. With the help of the transverse projector, A'^ = g^^ — u^u^ , 
we consider 

T'^ = euV + q'u^ + u'q^ - {p + H) A*^' + 7r*^ (1) 

where we define the transverse energy flow four vector (or heat flow in case one uses Eckart 's definition), as = 
A^T'^^Uj, the hydrostatic pressure, p, the bulk viscous pressure 11, where (p -|- 11) = — ■|AijT*-', and the shear stress 
tensor, tt*^ = T<'^> - {p + U)A'K The T<'^> = i(r'=" r"'=)Aj.A^ - iA'^ AfenT*^" notation stands for a particular, 
symmetrized and traceless combination of indices reflecting the same property of the shear stress tensor, TTy . 
The energy flux four vector is hence given by 

E^ = eu' + q\ (2) 

while its Lorentz invariant square length becomes 

E,E' = e^ + q,q\ (3) 

We note that since g^u' = and is a timelike vector, is spacelike, so g^q* < 0. The Lorentz scalar square of the 
energy-momentum tensor reads as 



(4) 
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where the number 3 = A* reflects the space dimensions. Our present ansatz generahzing the local invariant energy 
density including irreversible processes is given by: 

L = E,E' + {E,E' - Ty-T'^) + {E,E' - e^). (5) 

Here the first two terms, 2EiE'^ — TijT^^ — — 3(p + H)-^ — n^^TTij, contain the square of the energy density and the 
dissipative contributions due to the hydrostatic pressure plus bulk pressure and shear stress. The last term, EiE^ — e^, 
returns the square of the absolute value of the heat flow. Substituting these values our flducial scalar entering the 
equation of state is given by the square root of 

L = e^ +q,q' -3{p + n)2 - tt.jTt'^ . (6) 



We interpret this construction as the following physical picture: The effective energy density e = ^ L(p = 0) is equal 
to the familiar one, e, in the absence of dissipation. Dissipation decreases e compared to e, and since the entropy is a 
monotonic growing function of internal energy for systems with positive absolute temperature, the entropy is maximal 
at no dissipation. This construction is akin to the Miiller and Israel-Stewart approach in its spirit The p ~ 
version of eq. ([6]) separates the static from the dissipative parts of the pressure, according to the basic presumption 
of these theories. In our approach the equation of state is given by a particular function, 



s{e,n) = s |^^e2 + g^g'' - 3n2 - n^jn'o^nj = s (e, g', H, tt'-' , n) . (7) 

For thermodynamical systems without heat conduction and dissipative effects, = 0, 11 = and tt^ = 0, the 
traditional equilibrium relation emerges, s(e, 7i)|^i^Q_n=o,ir*-'=o = s(e,n) — s(e, 0', 0, O*-' , n). 

Moreover, s(e, n) achieves maximum at this point; so it is ensured that near equilibrium the dissipative currents 
relax. An expansion of e for small dissipative currents to leading order leads to an Israel-Stewart type of ansatz, 
however with fixed coefficients of the quadratic terms, 

s(6,n) « s(e,n) + 1 {q,q^ - - 7r.,7r'^) + ... (8) 

Here identifying the inverse equilibrium temperature via l/Teg = ds{e,n) / de , the Israel-Stewart coefficients are given 
by /3o = 3/e for the bulk viscosity, /3i = 1/e for the shear viscosity and (32 = 1/e for the heat fiow. This approach 
differs from the one obtained in the framework of kinetic theory, nevertheless for an ideal relativistic gas (e = 3p) the 
coefficients are all inversely proportional to the pressure, and our result comes close to some of the coefficients in Ref. 
p,9i] . Furthermore our formula can be applied for matter containing massless particles, while some results calculated 
from relativistic kinetic theory in Ref. [l9| diverges for m = 0, as well as for vanishing pressure. Here we note that one 
might introduce different scalar functions in front of all newly introduced scalar terms when constructing the effective 
energy density, for example in order to match the thermodynamical coefficients introduced by Israel and Stewart. In 
the present paper we do not introduce such coefficients in order to keep the investigations simple and transparent. 

In the following we study the class of relativistic equation of states involved in eq. ([7]) without, surpassing this way 
the usual second order fluid dynamics of Israel-Stewart. We note that it is useful to introduce a common notation 
for all dissipative modifications. We consider = — D^, where = ^ <Zi<Z* + SH^ + Tiijir^^ from which the 
approximations, e « e — D^/{2e) and s(e,n) « s{e,n) — D^/{2eT) holds. 

The Gibbs relation can be obtained by inspecting the total differential of the entropy density ([7]) . We obtain two 
sets of intensive variables differentiating the indirect function, s(e(e, g^, H, 7r,j), n). The differential with respect to e 
and n, 

, , ds ds 1 fi 

as(e, n) = —de + —dn — -de — —dn , (9) 
oe on 

defines the effective temperature 9 and effective chemical potential /t by simple expressions. The introduction of 
these new intensive variables naturally mimics the standard thermodynamic relations in case of equilibrium, the 
interpretation of these quantities becomes clear in the equilibrium limit. 
On the other hand, since e, depends on e, g*, H and tt*-' , 

ds{e,q\Ii,TT'\n) = ^(-de^^dq'-—dIl-^dnA+^dn (10) 
oe \e e e e J on 

^ Lde + ^dq' - — dn - ^d^'^" - ^dn , 
T eT ^ eT eT T ' 
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enables us to introduce intensive variables associated to de, dq^, dll, dn"^^ and dn, respectively. As we will see, 
these intensive quantities appear in the diffusion, pressure and heat conduction relaxation equations, driving the 
thermodynamic system toward equilibrium. Therefore T is called equilibrating temperature distinct from the effective 
temperature 9 appearing in the state functions. Using the above relations we can easily establish the following relations 
between those quantities 

eT — e6, and e^ = efl, (11) 

where the functions 9, T and fi, /i are respectively equal and reduce to the equilibrium temperature and chemical 
potential, T^q = T = 9 and fi^q = = Aj when the energy flux and the viscous pressure vanish (g* = 0, 11 = 
0, TT*^ = 0). The equilibrium entropy density corresponds to the previously introduced s{e,n) = Seq{e,n), because 

de = TeqdSeq + Heqdn. 



III. THERMODYNAMIC PRESSURE 



Since the derivative with respect to proper time in a frame comoving with the local flow is given by d/dr — Uid'', 
the above total derivative multiplied with the temperature can be expressed as 

ds 

T— = Tdi(su') - Tsd,u\ (12) 
dr 

This leads to the following entropy balance equation 

T^^{su') = ^^{eu') + ^d^iq^u') - —d,{Uu') ~ ^a,(7r^''=u') - ^a,(nu') + pd^u\ (13) 



Here each term expresses a contribution to the entropy production as a product of an intensive parameter and the 
divergence of an extensive one. In particular the last term in eq. (|13p contains the mechanical work on a changing 
volume (considering dV/dr — VdiU^). The coefficient p turns out to be 

p = Ts{e, n) + fin . (14) 

e 

An equivalent way reminding the familiar thermodynamic relation is given as 

—p = 9s{e,n) + fm-e. (15) 

Using the comprised notation for the sum of dissipative modifications in the energy density, = — , while 
expanding the entropy density around the traditional formula, s(e,n) — Seq — D"^ /{2eTeq), an approximate formula 
arises for the pressure 

P W TeqSeq + pieqn ~ ^ + ' (16) 

On the other hand for the - till this point unspecified - parameter p in the energy-momentum tensor in the dissipative 
case we suggest to use the standard expression, 

p — Ts + fin — e. (17) 

With this definition in general it follows that, p — p + j e « Peq + D"^ / {2e), with p^q — T^qS^q + fieqn — e being the 
isotropic equilibrium pressure in the absence of dissipation. 

At this point it is a delicate question what part of the total pressure will be assigned to dissipative and what 
to non-dissipative effects; the thermodynamic interpretation of the parameter p occurring in the energy-momentum 
tensor has to be determined. While in the absence of dissipation p coincides with the isotropic, equilibrium hydrostatic 
pressure, Peg, from the analysis of the mechanical work done on the expanding volume the quantity p is relevant. We 
motivate our choice by the following brief analysis of the example of pure radiation. 

Our starting point is the following equation of state for an ideal gas of massless particles (radiation): 



s(e) = s(e,g\n,7r^^) = ae^/\ 



(18) 



5 



where a = (3a/4) is the Stefan-Bohzmann constant. The inverse temperatures in this case are given by 

This means that the effective energy density is, e = cr6'^, while the pressure from eq. (|14p becomes 

le2 1 T 

From this at no dissipation the famihar equation of state p^q — e/3 arises. From the energy balance one knows, that 
the cooling of an expanding system is driven by the quantity h = e + p. For the pure radiation this becomes 

4 £2 

h = e+p^Ts^ = 4p. (21) 

3 e 

From this study we conclude that the interpretation of e as internal energy must be developed consequently. In our 
understanding the thermodynamic pressure is given by p in eq. (|14p . 

IV. ENTROPY PRODUCTION AND LINEAR RESPONSE 

In order to proceed further by determining the entropy production we need the energy-momentum balance, 

djT'^ = 0. (22) 

This way using eq. ([T]), the energy balance can be expressed as 

u.djT'^ = d^{eu') + d,q' + u,q' + {p + Ii)d^u' - ir'^d^Uj = 0. (23) 

Here the over-dot denotes the proper time derivative, x — dx/dr — Uid^x. Using the fact that g^M* = we replace the 
term Uiq^ = —qiii^. The Euler-equation is given 

A* t^djT^^ = eii' + q'djU^ + q^dju' + A* t,q^ - A* ^d^P^^ = 0, (24) 

where P*^ = {p + H) A'^' + ir^i . 

Substituting the energy balance into the entropy balance (|13p we obtain 

Tdi(su') = ~^^q' + qiii' - U^^u' + TT'^diUj - fidi(nu') + ^q' - — H - ^r'. (25) 

e e e 

The particle four current is defined generally as — nu^ + z/' , where is the particle flux, which is spacelike in 
the local rest frame. Associating the chemical potential to a conserved particle, S^TV* = di{nu^) + diV^ ~ 0, and it 
allows to use di{nu^) = —diV^ in the entropy balance. Therefore, substituting eq. (I^S)) into eq. (I^S)) we get 



where according to the Second Law of thermodynamics the entropy production is non-negative. We can see, that 
convenient definition of the entropy flux is 
Hence the entropy four current is defined as 











T 


~Y 






T 





s{e, n^u' + ^-'-^ + {q,q^ 7r,.^^'=) (27) 

Here the first term in the last line denotes the entropy of a perfect fluid carried by the flow, the second term denotes 
the entropy flux due to heat and particle flux. In case one uses Eckart's deflnition for the flow field, i/"^ = 0, the energy 
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flux (f equals to the heat flow. On the other hand using the definition of Landau and Lifshitz for the flow, (f = 0, the 
energy flux vanishes in the local rest frame and heat flux is defined by P = —{e + p)/nv^. In a baryon free matter 
one obviously has 1/^ = 0, therefore there is no temperature equilibration by heat conduction. 

The construction involved in (|27p truncated terms up to linear order in dissipative quantities corresponds to Eckart's 
first order theory of relativistic dissipative fluid dynamics. The introduction of quadratic terms are generally referred 
to as second order theories of dissipative fluid dynamics [lll.[20l|. Our approach (cf. first line in cq. ()27p ) contains an 
infinite series of higher order terms. 

A complete set of second order terms in the entropy four current of relativistic fluids was proposed by Israel and 
Stewart [1^, with coefficients /3o,/3i and (32, for the quadratic terms in 11, and tt'^ respectively. They introduced 
also OfQ Eind ai coefficients in front of viscous-heat flux coupling terms like Hq^ and Tr'-^g^. Since, our main ansatz is 
based on Lorentz scalar second order quantities involved in the dissipative relativistic energy density, it is clear that 
such vector terms do not appear in our present approach. 

The Miiller, Israel and Stewart method modifies and generalizes explicitly the entropy four current of Eckart, 
through which the definition of the local rest frame entropy density is extended. They give the most general isotropic 
and second order expression in a pure mathematical way. Our method is based on constructing local corrections to 
the energy density, and then define the entropy four current along the lines as done by Eckart. However, expanding 
the corrections to the equilibrium entropy density (|27p results in a formula resembling the definition of Israel and 
Stewart. 

The nonnegativity of the entropy production in (|26p can be ensured if each term is separately non-negative. There- 
fore in case of isotropic materials the coefficients of the energy flux q^ (orthogonal to u', i.e. the heat conduction), the 
bulk and shear viscosity terms multiplying 11 and TTy in the above expression and the chemical diffusion contribution 
are treated as being proportional to the corresponding coefficient in the linear response approximation. This leads to 
the following equations: 

q' = -A (tA'Hj + ^^"^j + V'T^ , (28) 

n = -c (d,u' + ^n) , (29) 

TTy = 277 (^d<iUj> - ' (^0) 

z.^ = -<yd'^, (31) 

where the short hand notation V* = A'^9j has been introduced. The heat conductivity A, the bulk viscosity C 
and the shear viscosity rj are non- negative transport coefficients. The diffusion term (last line in eq. (|3ip ) and the 
heat conduction term (first line) may in general show cross-couplings (cf. Soret effect t21j). With the above simple 
construction, we arrive at the following formula for the entropy production, 

In the following discussions we concentrate on the heat conduction and viscosity evolution and neglect particle 
diffusion phenomena. Based on the above considerations we obtain the following evolution equations for the viscosity 
terms: 

^A^^q^ + lw^T + A^^u. + j^q^ = 0, 

-il + d.u' + lu = 0, 
e C 

-TTy - d^iUjy + — TTy 0. (33) 

The corresponding relaxation times are hence given by Tq = XT(3i = XT/e for the heat conduction, by Tti = /?oC = 3C/e 
for the bulk viscosity and by — 277/32 = 2r]/e for the shear viscosity. These relaxation times are perfectly finite even 
for massless matter. 

The above relaxation equations closely resemble the truncated Israel-Stewart form. Such forms of the relaxation 
equations are extensively utilized in quark matter research 0, [2^ . [23j . 

The main difference to our approach lies in terms which contain the space-time derivative of the thermodynamic 
coefficients multiplied by the flow. As argued by Heinz et. al. [l^l, these Israel-Stewart quantities are rather small. 
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and one can neglect them safely. In the case they were large, they become unbounded and the system may destabilize. 
In our approach such terms, containing derivatives of the Israel-Stewart coefficients, do not appear. However, these 
equations already capture the essential features of relaxation phenomena and we expect that it will result in a causal 
and stable theory. 

These relaxation equations can be solved parallel to the energy flux equation (j23|) and the Euler equation describing 
the evolution of the flow u*. 



V. GENERALIZED BJORKEN FLOW 



grad = N,d' = 7 ( + ) . (35) 



Let the basis of coordinates be given by the tetrad ejj, such that = tcq + re^, and therefore 
dx^ = dreg + rdrjCi + dre'2 + rd(t>e\ . The unit four vectors satisfy the orthogonality relations 

■i j 

eleig^' = 9'' (34) 

with the Minkowski metric tensor gij = diag(l, —1,-1,-1). The flow velocity field and an orthogonal spacelike field 
are given by 

= 7(4 + «ei) > 

iV' = 7 {vel + el) , (35) 

where 7 = l/Vl — w^, UiU^ — +1, NiN^ = —1 and u^u^ — N^N^ = CqCq — e\e\ so m* and A^* span the same two- 
dimensional subspace of the spacetime as Cg and e\ . It is therefore natural to consider derivatives in the direction of 
these vectors, we call them as 'dot' and 'grad': 

^, f d V d 

\OT T orj 
d_ Ic 
dr T di] 

Four-divergences of the flow and its orthogonal are given as 

T 

B = a^u* = - +72grad(u). (37) 

T 

All partial derivatives can be expressed by these quantities: 

d'u^ = {Au' ~ BN') N^, 

d'N^ = {Au'-BN')u^. (38) 
The transverse projection tensor in our basis is given by 

A'^ = g'^ - u'u' = - (^N'N^ + e^ei + ele(^ . (39) 

It has the properties UfcA*^' = and A^fcA*^* — N^. This helps to obtain the traceless symmetric part of the derivative 
tensor of the velocity field, important to entangle the shear pressure term: 

= B (^-N^N^ - Ia^^^ - B (^^^N^N^ + y^ei + ie?,e|) . (40) 

The shear pressure tensor tt*-' is proportional to the expression in the bracket, we use a factor 3/2 in the definition 

= Ar, V) (-N'N= + ^e^ei + ^e^^] . (41) 
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Now n'^^TTij = 37r^/2 follows. The usage of the {u\ A^*) reference frame allows us to give the energy momentum tensor 
^ as follows: 

q{N'u^ + u'N^) + aN'N^ + (3 {ele{ + e*e^) , (42) 

with a = p + n — TT and (3 = p + 11 + 7r/2. It consists of projector terms, but the term proportional to g, which is due 
to the — q{T, ri)N^ form. 

We consider the energy momentum conservation, including all above terms in the diT^^ — equation. Its general 
form is given by 

u,dot(r^) = Ar,grad(r^), (43) 

so all equations describe a balance between 'dot' and 'grad' terms. In the case of the Bjorken flow these operations 
coincide with the time and rapidity derivatives, but for a longitudinal accelerating flow ansatz not. Together with 
the relaxation equations due to the linear response assumption we obtain the following set of dynamical equations 
corresponding to (USD, (EH), dM]), (EH) and ^ respectively: 

dot(e) + (e + a)B + 2qA + grad((2') = 0, 

dot((7) + (e + a)^ + 2(7B + grad(a) = 0, 
AT 

— dot((?) + q + XT A + Agrad(T) = 0, 
e 

— dot(n) + n + cs = 0, 

e 

— dot(7r) + TT - = 0. (44) 

Here the first equation describes the cooling due to expansion, the second is the Euler equation describing the 
acceleration of the flow due to pressure gradients, the third is the Fourier heat conduction equation supplemented 
with a relaxation term, while the fourth and fifth equations describe the relaxation of bulk and shear viscosity. 

It is interesting to consider that class of solutions when the quantities under investigation depend only on the time 
variable, r. Then denoting hy f — df/dr for such functions wc arrive at 

e = -{e + a)B - 2qA - vq, 
q = —{e + a)A — 2qB — va, 

- -AT^'^-^^-T^' 

e e ~ 

with 

V ~ 1 

A = - + 7^{), and B = - + -f'^vv. (46) 

T T 

In this form the cooling and the Eulcr equation show a quite symmetric role in the evolution, and furthermore it seems 
that one would obtain two equations for q. However, the Euler equation has to be used to describe the change of the 
flow, so it has to be regarded as an equation for v included in the variables A and B. The u = assumption reveals 
that A — ii and B — 1/r, so indeed the roles are not symmetric. In this case we obtain from the Euler equation 



T ^ 

and then using the Fourier equation for eliminating q from the above result we arrive at 



2 



a \A7 T 



1 . 



(48) 



One concludes that only for g = can the v = Q Bjorken flow remain stationary. This on the other hand consequently 
solves the Fourier equation. Hence we pointed out, that the energy flux, can be consequently set to zero when 
considering the stationary Bjorken flow independent of the Eckart or Landau-Lifshitz choice of the flow frame. 
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VI. EXPANSION, COOLING AND RE-HEATING IN THE SCALING SOLUTION 



It is customary to investigate the relaxation of viscosity by utilizing a stationary solution of the non-dissipative 
system, in particular for the quark gluon plasma evolution the Bjorken flow pattern [2^. We have seen already that 
the only assumption consequent with a stationary Bjorken flow (v = 0) is q = 0. In this case the acceleration of the 
flow and the heat current remain zero, due to eqs. (|47l [48]) . while the energy density, the bulk and shear viscosity 
relax in a coupled manner. Here, eqs. simplifies to 

e + p TT — n 



^ e e 

e 2e 

TT H TT = . 

2r] 3t 



(49) 



The corresponding results of the first order (Navier-Stokes) hydrodynamics are obtained by neglecting the tt and 11 
terms in the above equations. Then one considers tt^^^ — Arj/{?>T) and n*^^^ — —C,/t and observes a re-heating of the 
expanding quark gluon plasma: 

4 477 C 

In particular for early enough times the terms on the right hand side, being proportional to l/r^, dominate the 
evolution [13, [111- 

For the radiative equation of state, s = Icr^/'^e^/^, and e + p 4e^/(3e). According to (gl]), = - 311^ - Zt:'^ /2, 
thus e+ p — 4e/3 — 4(11^ + 7r^/2)/e. Therefore, the cooling due to the expansion is reduced, the dissipative terms 
physically re- heat the system. There is nothing artificial about it. This is fact more apparent for the proper time 
derivative of the effective internal energy density together with eqs. (j49[) . where the terms Il/r and tt/t do not appear 
explicitly, 

3t e V 47; C / 

This form of the cooling equation reflects the fact that the source of physical re-heating consists of the dissipative 
terms only, quadratic in dissipating momentum fluxes and inversely proportional to the linear response coefficients. 
For the sake of simplicity in the following we neglect the bulk viscosity, therefore eqs. ([^5]) reduce to 



4 37r2 / 3 /7r\2 



These equations augmented with the radiation EOS are to be solved. In the following, we will show the numerical 
results for the above equations in case of an ideal QGP with 3 massless quarks and 16 gluonic degrees of freedom, where 
the Stefan-Boltzmann coefficient is given as = 47.57r^/30, while the coefficient of viscosity is t; = 7/os = 7/o|(o'e^)^/*. 
Without viscous pressure terms, for and equilibrium EOS, one has to replace e with e. The initial conditions are given 
following Ref. ^27j. The starting time for the hydrodynamical evolution is, tq = 0.6 fm/c, the initial energy density 
eo = Co = 30 GeV/fm^, the initial viscosity to entropy density is 7;o = 0.4, while the initial shear for , ttq = 0. Here 
we note that other initial conditions and initial values are also possible, see for example f28|,[29|. 

To calculate the temperature eqs. ([5 ^ - ([55)) are solved together with the standard QGP EOS in four special cases. 
We give solutions for a perfect fluid, where tt = 0; for a first order dissipative fluid, i.e., the relativistic Navier-Stokes 
equations (NS), where tt = 477/3t and e = e; for second order dissipative fluid, i.e., Israel-Stewart type fluid (IS), 
where we set e = e and four our higher order equations (HO) without any of the previous simplifications. The initial 
shear for the transport equation is, ttq = 0. The corresponding thermodynamic temperatures are denoted by Tiu for 
a perfect fluid, T^rg for the NS equations, T/5 for the IS type equations and T^q and Oho for the two temperatures 
of our higher order theory. The existence of these temperatures is the outcome of the thermodynamical classification 
and treatment as given in eqs. ^ [TO]). Their different physical role is clear from the structure of the theory; T is 
responsible for the equilibration and appears in the EOS. 
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0.15 T [fm/c] 

2 4 6 8 10 

FIG. 1: The evolution of temperature as function of the proper time for the initial conditions and equations given in the text. 
The dotted line shows the cooling of a perfect fluid, Tjd, the thin dashed and thin full lines are for the NS and IS equations, 
denoted by Tns and Tjs- The thick dashed and thick full lines are for the two temperatures of the HO equation, respectively 
denoted by Tho and Oho- 

On Fig. [U the evolution of the temperature is shown for the different simplified models. The overall behavior of 
the results is closely the same. For a longer time cooling is fastest for the ideal fluid, without dissipation, slower 
for the NS and IS fluids and slowest for the HO fluid, due to the reduced value of the expansion strength, h — 
4e/3 — 4(n^ + 7r^/2)/e, compared to the case when, h — 4e/3. The equilibrating temperature T and the effective 
temperature 9 tend to each other with the decreasing dissipation. For the HO solution, the equilibrating temperature 
decreases faster since T = Oe{e^ + 1.57r2)-i/2 and tt initially rises then decreases to the Navier-Stokes limit. 

To phenomenologically study the stability of the system, investigation related to the Reynolds number is standard 
practice. The inverse Reynolds number, = 7r/(e + p), is the ratio between dissipative and non-dissipative 

quantities, i.e., the ratio of dissipation to the strength of expansion. The equation for the energy density in (|49p can 
be re-written in the following form, 

e=-(l-i?-^), (54) 



T 



where h = e + p. From eq. (|5ip we get 



where 



e = - (1 - R-') , (55) 



T 



= (56) 

We easily realize that in both cases the energy increases as long as the inverse of the Reynolds number is smaller than 
one. This leads to phenomenological upper bounds for the dissipative pressure, i.e., tt < Ap. 

For the first order theory of Eckart one can show that in case the dissipative pressure becomes lar ger, init ially 
or otherwise, than four times of the isotropic pressure, the solution to the equations becomes unstable |20l. [3ll. |32|| . 
There reheating is closely related to stability. According to our knowledge the upper bound on dissipative quantities 
is not explicitly or quantitatively well specified for higher order theories, therefore their domain of applicability and 
the stability conditions of the Bjorken flow and its generalizations are not clear. 

On Fig. [21 the evolution of the inverse Reynolds number is shown for the two definitions given before. The 
dashed line is for the NS equations while the full lines for the IS equations. In the NS case the inverse Reynolds 
number decreases from its initial value and asymptotically approaches the perfect fluid limit due to the decrease of 
the expansion rate. For the IS and HO solutions, the system initially cannot compete with the fast expansion and 
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first departs from equilibrium, and only later relaxes to the NS solutions. Finally, we may conclude that there is no 
reheating in any of above presented cases because the inverse of the Reynolds number is less than one. 

Because our modified Reynolds number is explicitly time dependent, the separation of cooling and reheating solu- 
tions is not straightforward without the solutions of the equations. In fig. [3] we compare the reheating capabilities 
of the different theories. The minimal initial energy density for a reheating solution is plotted as function of starting 
proper time in case of zero initial viscous shear pressure. The shear viscosity was 770 = 0.3 in the left figure and 

rjo — 0.08 at right. For the first order theory of Eckart one can give the explicit condition as eo — cr (^j^^ ■ For 

the Israel-Stewart and four our higher order theory the corresponding curves can be well approximated as Cq — bT^'^, 
where the b parameter values are tabulated in table 1. 



Va 


Eckart 


Israel- Stewart 


This paper 


0.3 


6 • lO-'' 


5.6-10-7 


2.67- 10-4 


0.08 


3 • 10-6 


2.89- 10-9 


1.75-10-4 



TAB. 1: b [GeV fm/c^] parameter values characterizing the conditions of reheating for the different dissipative fluids. 



eo[GeV/fm^] 
10 




eo[GeV/fm^] 
10 



O.I 0.2 0.5 I.O 2.0 



To[fm/c] 




0.5 1.0 2.0 



To[fm/c] 



FIG. 3: The initial energy density to get a reheating flow as a function of the starting proper time for rjo = 0.3 to the left 
and ?7o = 0.08 to the right. The thin dashed and thin full lines are for the NS and IS fluids and the thick full line is for the 
HO model. Initial conditions below the lines lead to reheating solutions. The dots indicate initial conditions considered to be 
reahstic for the RHIC and LHC experiments [27l.l28|. 
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VII. DISCUSSION 

In our approach we have fixed the thermodynamic coefRcients, akin to the Pi coefficients of the Israel-Stewart 
theory, respectively or the a coefficient of the Ottinger-Grmela theory [H, [sj to particular values, that otherwise 
could have been calculated from kinetic models. There are several reasons why a direct phenomenological approach 
can be fruitful in heavy ion physics. 

1. The available transport coefficients calculated from kinetic theory are not allways realistic and sometimes are 
controversial. Not realistic, when they are related to oversimplified microscopic properties (like one component 
ideal gases) fioj. In this approach some structural properties are considered for the second order coefficients 
[23]. 

2. Certain results of hydrodynamic calculations are independent of the exact values of the second order coefficients: 
A difference in the initial values of energy density, equilibrium and dissipative pressures can be compensated 
by an appropriate choice of the initial proper time (tq) in the Bjorken flow scenario. We have argued that first 
order theories cannot be excluded by causality reasons, because the actual violation of causality may be beyond 
the validity range of the hydrodynamic approach [l3|. Baier et. al. has been arguing that the difference of 
the solutions of a first order and a second order theory do appear at microscopic distances, beyond the validity 
range of hydrodynamics 0]. Several theories of relativistic dissipative fluids (e.g. Israel-Stewart, the divergence 
type Geroch or Miiller-Ruggieri theories, Baier et al., Koide et al., Ottinger-Grmela) may belong to this class. 

However, it is easy to see, that generic instabilities of the corresponding theory could destroy the above argu- 
mentation. Up to know the instability of the first order theories and stability conditions in Israel-Stewart theory 
[Tol . [30| has been well investigated and for the 7r = n = 0, g = minimal version of our higher order approach 
were pointed out p^ . According to our knowledge, the conditions of stability of homogeneous solutions are not 
known for the rest of the above mentioned theories. 

3. There is no a priori reason why one should end at second order extension of the first order theory, beyond 
convenience. 

Within our thermodynamic framework, not showing generic instabilities for ^ we have investigated the role of 
heat conduction together with a generalized, longitudinally accelerating Bjorken fiow in more detail. Our conclusion 
was that for stationary flow the heat flux must vanish also in a general frame (also in the Eckart frame), there is no 
need for the the customary Landau-Lifshitz condition (g* = 0). On the other hand the reheating effect is a consequence 
of dissipative physical phenomena, not of generic instabilities. 

A peculiar property of our approach to the relativistic internal energy is the distinctness of the derivatives of the 
entropy with respect the total or internal energies. This fact resulted in a doubled set of non-equilibrium intensive 
thermodynamic quantities according to the respective derivatives (T and 9, fi and //, p and p). We have seen that 
the intensive variables related to the internal energy (6,(i,p) resulted in more natural thermodynamic relations: The 
intensive variables related to the total energy {T,ii,p and the coefficients of the viscous pressure terms in the Gibbs 
relation (jlOp ) are rather related to the energy equilibration due to dissipative processes. In particular the gradient of 
T appeared in the generalized Fourier law (|28p but for the radiation thermodynamics it was most straightforward to 
consider 9 as temperature. 
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